1 Introduction
This webbook contains all the code used for data analysis in study of the individual-level metagenomic data of Podarcis muralis and Podarcis liolepis lizards from different environments during an experimental setup.
1.1 Prepare the R environment
1.1.1 Environment
To reproduce all the analyses locally, clone this repository in your computer using:
RStudio > New Project > Version Control > Git
And indicating the following git repository:
Once the R project has been created, follow the instructions and code chunks shown in this webbook.
1.1.2 Libraries
The following R packages are required for the data analysis.
# Base
library(R.utils)
library(knitr)
library(tidyverse)
library(devtools)
library(tinytable)
library(janitor)
# For tree handling
library(ape)
library(phyloseq)
library(phytools)
# For plotting
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(ggnewscale)
library(gridExtra)
library(ggtreeExtra)
library(ggtree)
library(ggh4x)
# For statistics
library(spaa)
library(vegan)
library(Rtsne)
library(geiger)
library(hilldiv2)
library(distillR)
library(broom.mixed)
library(corrplot)
library(nlme)
library(pairwiseAdonis)2 Prepare data
2.1 Load data
Load the original data files outputted by the bioinformatic pipeline.
2.1.2 Read counts
read_counts_raw <- read_tsv("data/genome.count.tsv") %>%
rename(genome=1)
#Transformation of read_counts to combine data from both sequence rounds
merge_and_rename <- function(read_counts_raw) {
read_counts_raw %>%
# Gather the columns into long format
pivot_longer(cols = -genome, names_to = "col") %>%
# Extract prefix
mutate(prefix = gsub("^(.*?)_.*", "\\1", col)) %>%
# Group by prefix and genome, then summarize
group_by(prefix, genome) %>%
summarize(value = sum(value)) %>%
# Spread the data back to wide format
pivot_wider(names_from = prefix, values_from = value)
}
read_counts <- merge_and_rename(read_counts_raw)2.1.3 Genome base hits
genome_hits_raw <- read_tsv("data/genome.covered_bases.tsv") %>%
rename(genome=1)
#Transformation of genome_hits to combine data from both sequence rounds
merge_and_rename <- function(genome_hits_raw) {
genome_hits_raw %>%
# Gather the columns into long format
pivot_longer(cols = -genome, names_to = "col") %>%
# Extract prefix
mutate(prefix = gsub("^(.*?)_.*", "\\1", col)) %>%
# Group by prefix and genome, then summarize
group_by(prefix, genome) %>%
summarize(value = sum(value)) %>%
# Spread the data back to wide format
pivot_wider(names_from = prefix, values_from = value)
}
genome_hits <- merge_and_rename(genome_hits_raw)2.1.5 Genome quality
genome_quality <- read_tsv("data/quality_report.tsv") %>%
select(
genome = Name,
completeness = Completeness,
contamination = Contamination,
length = Genome_Size,
gc = GC_Content
)
genome_quality<-genome_quality %>%
mutate (genome=str_remove_all(genome,".fa"))
#Filter MAGs with over 70% completeness and less than 10% contamination
genome_quality <- genome_quality %>%
filter(completeness > 70 & contamination < 10)2.1.6 Genome tree
genome_tree <- read_tree("data/gtdbtk.backbone.bac120.classify.tree")
genome_tree$tip.label <- str_replace_all(genome_tree$tip.label,"'", "") #remove single quotes in MAG names
#Filter genome_taxonomy to keep MAGs with over 70% completeness and less than 10% contamination
genome_taxonomy <- genome_taxonomy %>%
semi_join(genome_quality, by = "genome")
genome_tree <- keep.tip(genome_tree, tip=genome_taxonomy$genome) # keep only MAG tips2.2 Create working objects
Transform the original data files into working objects for downstream analyses.
2.3 Filtering
4 samples are removed from the analysis due to their low sequencing depth.
#Counts_raw
columns_to_exclude <- c("AD16","AD23", "AD25", "AD91") # Columns to exclude
read_counts <- read_counts %>%
select(-columns_to_exclude)
#Metadata
sample_metadata <- sample_metadata %>%
filter(Tube_code != "AD16") %>%
filter(Tube_code != "AD23") %>%
filter(Tube_code != "AD25") %>%
filter(Tube_code != "AD91")
#Coverage_table
columns_to_exclude <- c("AD16", "AD23", "AD25", "AD91") # Columns to exclude
genome_coverage <- genome_coverage %>%
select(-columns_to_exclude)##Removal of samples of individual LI1_2nd_6
#Counts_raw
columns_to_exclude <- c("AFV13","AD06", "AD31", "AE03", "AF12") # Columns to exclude
read_counts <- read_counts %>%
select(-columns_to_exclude)
#Metadata
sample_metadata <- sample_metadata %>%
filter(Tube_code != "AFV13") %>%
filter(Tube_code != "AD06") %>%
filter(Tube_code != "AD31") %>%
filter(Tube_code != "AE03") %>%
filter(Tube_code != "AF12")
#Coverage_table
columns_to_exclude <- c("AFV13","AD06", "AD31", "AE03", "AF12") # Columns to exclude
genome_coverage <- genome_coverage %>%
select(-columns_to_exclude)2.3.1 Filter reads by coverage
#Filter read_counts for the MAGs with over 70% completeness and less than 10% contamination
read_counts <- read_counts %>%
semi_join(genome_quality, by = "genome")
min_coverage=0.3
read_counts_filt <- genome_coverage %>%
mutate(across(where(is.numeric), ~ ifelse(. > min_coverage, 1, 0))) %>%
mutate(across(-1, ~ . * read_counts[[cur_column()]])) 2.4 Load data statistics
2.4.1 Raw reads
raw_reads <-
"data/multiqc_general_stats_2.txt" %>%
read_tsv() %>%
select(
sample = Sample,
raw_reads = `total_sequences`
) %>%
mutate(sample_prefix = str_extract(sample, "^[^_]+")) %>%
group_by(sample_prefix) %>%
summarize(raw_reads = sum(raw_reads, na.rm = TRUE)) %>%
rename(sample = sample_prefix) %>%
mutate(sample = str_remove(sample, "^fastp \\|\\s*"))2.4.2 Quality-filtered reads
fastp_reads <-
"data/multiqc_general_stats.txt" %>%
read_tsv() %>%
select(
sample = Sample,
trimmed_reads = `total_sequences`
) %>%
mutate(sample_prefix = str_extract(sample, "^[^_]+")) %>%
group_by(sample_prefix) %>%
summarize(trimmed_reads = sum(trimmed_reads, na.rm = TRUE)) %>%
rename(sample = sample_prefix) %>%
filter(!str_detect(sample, "nonlizard \\|")) %>%
filter(!str_detect(sample, "lizard \\|")) %>%
filter(!str_detect(sample, "refseq500 \\|")) %>%
mutate(sample = str_remove(sample, "^fastp \\|\\s*"))2.4.3 Host-mapped reads
host_mapped <-
"data/multiqc_general_stats.txt" %>%
read_tsv() %>%
filter(str_detect(Sample, "lizard")) %>%
select(
sample = Sample,
host_mapped = `reads_mapped`,
mapping_total = `raw_total_sequences`
) %>%
mutate(
host_unmapped = mapping_total - host_mapped
) %>%
filter(!is.na(host_mapped)) %>%
separate(
col = sample,
into = c("host_name", "sample"),
sep = " \\| "
) %>%
rename(mapped = host_mapped, unmapped = host_unmapped) %>%
select(-mapping_total) %>%
pivot_longer(-host_name:-sample) %>%
mutate(
name = str_glue("{name}_{host_name}")
) %>%
select(-host_name) %>%
pivot_wider() %>%
mutate(sample = str_extract(sample, "^[^_]+")) %>%
group_by(sample) %>%
summarize(
mapped_lizard = sum(mapped_lizard),
unmapped_lizard = sum(unmapped_lizard)
) 2.4.4 Prokaryotic fraction
singlem <-
"data/singleM.tsv" %>%
read_tsv() %>%
distinct() %>%
mutate(
sample = sample %>% str_remove_all("_1$"),
read_fraction = read_fraction %>% str_remove("%") %>% as.numeric(),
read_fraction = read_fraction / 100
) %>%
select(
sample,
singlem_prokaryotic_bases = bacterial_archaeal_bases,
singlem_metagenome_size = metagenome_size,
singlem_read_fraction = read_fraction,
) %>%
mutate(sample_prefix = str_extract(sample, "^[^_]+")) %>%
group_by(sample_prefix) %>%
summarize(
singlem_prokaryotic_bases = sum(singlem_prokaryotic_bases),
singlem_metagenome_size = sum(singlem_metagenome_size),
singlem_read_fraction = mean(singlem_read_fraction)
) %>%
rename(sample = sample_prefix)2.4.5 MAG-mapped reads
mag_mapping <-
read_tsv("data/contig.count.tsv.xz") %>%
pivot_longer(-sequence_id) %>%
summarise(value = sum(value), .by = "name") %>%
rename(sample = name, mapped_mags = value) %>%
mutate(sample_prefix = str_extract(sample, "^[^_]+")) %>%
group_by(sample_prefix) %>%
summarize(
mapped_mags = sum(mapped_mags)
) %>%
rename(sample = sample_prefix)2.5 Prepare color scheme
AlberdiLab projects use unified color schemes developed for the Earth Hologenome Initiative, to facilitate figure interpretation.
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
arrange(match(genome, genome_tree$tip.label)) %>%
select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
pull(colors, name=phylum)3 Data statistics
3.2 DNA fractions
#Overall
data_stats %>%
mutate(mapped_perc=mapped_mags/trimmed_reads) %>%
summarise(mean=mean(mapped_perc),sd=sd(mapped_perc))# A tibble: 1 × 2
mean sd
<dbl> <dbl>
1 0.438 0.207
data_stats %>%
mutate(
low_quality = raw_reads - trimmed_reads,
unmapped_reads = trimmed_reads - mapped_lizard - mapped_mags
) %>%
select(sample, low_quality, mapped_lizard, mapped_mags, unmapped_reads) %>%
pivot_longer(-sample) %>%
mutate(name=factor(name,levels=c("low_quality","mapped_lizard","unmapped_reads","mapped_mags"))) %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(!is.na(time_point)) %>%
ggplot(aes(x = sample, y = value, fill = name)) +
geom_bar(stat = "identity", position = "fill") +
scale_fill_manual(name="Sequence type",
breaks=c("low_quality","mapped_lizard","unmapped_reads","mapped_mags"),
labels=c("Low quality","Mapped to host","Unmapped","Mapped to MAGs"),
values=c("#CCCCCC", "#bcdee1", "#d8b8a3","#93655c"))+
facet_grid(~time_point, scales = "free", labeller = labeller(time_point=c("0_Wild"="Wild", "1_Acclimation"="Acclimation", "2_Antibiotics"="Antibiotics",
"3_Transplant1"="Transplant1", "4_Transplant2"="Transplant2", "5_Post-FMT1"="Post1",
"6_Post-FMT2"="Post2"))) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=0)) +
labs(y="DNA sequence fraction",x="Samples")3.3 Recovered microbial fraction
data_stats %>%
mutate(
unmapped_reads = trimmed_reads - mapped_lizard - mapped_mags,
mag_proportion = mapped_mags / (mapped_mags + unmapped_reads),
singlem_read_fraction = singlem_read_fraction
) %>%
select(sample, mag_proportion, singlem_read_fraction) %>%
mutate(
mag_proportion = if_else(singlem_read_fraction == 0, 0, mag_proportion),
singlem_read_fraction = if_else(singlem_read_fraction == 0, NA, singlem_read_fraction),
singlem_read_fraction = if_else(singlem_read_fraction < mag_proportion, NA, singlem_read_fraction),
singlem_read_fraction = if_else(singlem_read_fraction > 1, 1, singlem_read_fraction)
) %>%
pivot_longer(-sample, names_to = "proportion", values_to = "value") %>%
mutate(
proportion = factor(
proportion,
levels = c("mag_proportion", "singlem_read_fraction")
)
) %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(!is.na(time_point)) %>%
ggplot(aes(x = sample, y = value, color = proportion)) +
geom_line(aes(group = sample), color = "#f8a538") +
geom_point() +
scale_color_manual(name="Proportion",
breaks=c("mag_proportion","singlem_read_fraction"),
labels=c("Recovered","Estimated"),
values=c("#52e1e8", "#876b53"))+
theme_minimal() +
facet_grid(~time_point, scales = "free", space="free", labeller = labeller(time_point=c("0_Wild"="Wild", "1_Acclimation"="Acclimation", "2_Antibiotics"="Antibiotics",
"3_Transplant1"="Transplant1", "4_Transplant2"="Transplant2", "5_Post-FMT1"="Post1",
"6_Post-FMT2"="Post2")))+
labs(y = "Samples", x = "Prokaryotic fraction") +
scale_y_continuous(limits = c(0, 1)) +
theme(
axis.text.y = element_text(size = 4),
axis.text.x = element_text( angle = 90, vjust = 0.5, hjust = 1, size = 0),
legend.position = "right"
)3.3.1 Domain-adjusted mapping rate (DAMR)
data_stats %>%
mutate(
unmapped_reads = trimmed_reads - mapped_lizard - mapped_mags,
mag_proportion = mapped_mags / (mapped_mags + unmapped_reads),
singlem_read_fraction = singlem_read_fraction
) %>%
mutate(damr=pmin(1, mag_proportion/singlem_read_fraction)) %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(!is.na(time_point)) %>%
select(sample,damr, time_point, species, type) %>%
#group_by(type) %>%
#summarise(mean=mean(damr),sd=sd(damr)) %>%
tt()| sample | damr | time_point | species | type |
|---|---|---|---|---|
| AC79 | 1.0000000 | 1_Acclimation | Podarcis_muralis | Control |
| AC80 | 1.0000000 | 1_Acclimation | Podarcis_muralis | Treatment |
| AC81 | 0.9813013 | 1_Acclimation | Podarcis_muralis | Treatment |
| AC82 | 0.7490771 | 1_Acclimation | Podarcis_muralis | Treatment |
| AC83 | 1.0000000 | 1_Acclimation | Podarcis_muralis | Control |
| AC84 | 1.0000000 | 1_Acclimation | Podarcis_muralis | Control |
| AC85 | 0.8363777 | 1_Acclimation | Podarcis_muralis | Control |
| AC86 | 0.8800275 | 1_Acclimation | Podarcis_muralis | Control |
| AC87 | 0.9672854 | 1_Acclimation | Podarcis_muralis | Treatment |
| AC88 | 0.3142108 | 1_Acclimation | Podarcis_muralis | Control |
| AC89 | 0.9338962 | 1_Acclimation | Podarcis_muralis | Treatment |
| AC90 | 1.0000000 | 1_Acclimation | Podarcis_muralis | Control |
| AC91 | 1.0000000 | 1_Acclimation | Podarcis_muralis | Control |
| AC92 | 1.0000000 | 1_Acclimation | Podarcis_muralis | Control |
| AC93 | 1.0000000 | 1_Acclimation | Podarcis_muralis | Treatment |
| AC94 | 0.9159348 | 1_Acclimation | Podarcis_muralis | Treatment |
| AC95 | 0.9239904 | 1_Acclimation | Podarcis_muralis | Treatment |
| AC96 | 0.9425962 | 1_Acclimation | Podarcis_muralis | Treatment |
| AC97 | 0.9394886 | 1_Acclimation | Podarcis_liolepis | Hot_control |
| AC98 | 1.0000000 | 1_Acclimation | Podarcis_liolepis | Hot_control |
| AC99 | 1.0000000 | 1_Acclimation | Podarcis_liolepis | Hot_control |
| AD01 | 1.0000000 | 1_Acclimation | Podarcis_liolepis | Hot_control |
| AD02 | 1.0000000 | 1_Acclimation | Podarcis_liolepis | Hot_control |
| AD03 | 1.0000000 | 1_Acclimation | Podarcis_liolepis | Hot_control |
| AD04 | 0.8896057 | 1_Acclimation | Podarcis_liolepis | Hot_control |
| AD05 | 1.0000000 | 1_Acclimation | Podarcis_liolepis | Hot_control |
| AD07 | 0.9515379 | 1_Acclimation | Podarcis_liolepis | Hot_control |
| AD08 | 0.7045003 | 2_Antibiotics | Podarcis_muralis | Control |
| AD09 | 1.0000000 | 2_Antibiotics | Podarcis_muralis | Control |
| AD10 | 0.9598720 | 2_Antibiotics | Podarcis_muralis | Control |
| AD11 | 0.6230176 | 2_Antibiotics | Podarcis_muralis | Control |
| AD12 | 0.7487463 | 2_Antibiotics | Podarcis_muralis | Control |
| AD13 | 1.0000000 | 2_Antibiotics | Podarcis_muralis | Control |
| AD14 | 0.4536529 | 2_Antibiotics | Podarcis_muralis | Control |
| AD15 | 1.0000000 | 2_Antibiotics | Podarcis_muralis | Control |
| AD17 | 0.8774871 | 2_Antibiotics | Podarcis_muralis | Treatment |
| AD18 | 0.4991585 | 2_Antibiotics | Podarcis_muralis | Treatment |
| AD19 | 0.4060523 | 2_Antibiotics | Podarcis_muralis | Treatment |
| AD20 | 0.5227482 | 2_Antibiotics | Podarcis_muralis | Treatment |
| AD21 | 0.6551882 | 2_Antibiotics | Podarcis_muralis | Treatment |
| AD22 | 0.9106725 | 2_Antibiotics | Podarcis_muralis | Treatment |
| AD24 | 0.5972915 | 2_Antibiotics | Podarcis_muralis | Treatment |
| AD26 | 0.9204825 | 2_Antibiotics | Podarcis_liolepis | Hot_control |
| AD27 | 0.5215059 | 2_Antibiotics | Podarcis_liolepis | Hot_control |
| AD28 | 0.8180607 | 2_Antibiotics | Podarcis_liolepis | Hot_control |
| AD29 | 0.7175978 | 2_Antibiotics | Podarcis_liolepis | Hot_control |
| AD30 | 0.9937574 | 2_Antibiotics | Podarcis_liolepis | Hot_control |
| AD32 | 0.6710817 | 2_Antibiotics | Podarcis_liolepis | Hot_control |
| AD33 | 0.6778368 | 2_Antibiotics | Podarcis_liolepis | Hot_control |
| AD34 | 0.4642326 | 2_Antibiotics | Podarcis_liolepis | Hot_control |
| AD36 | 0.8288365 | 3_Transplant1 | Podarcis_muralis | Control |
| AD37 | 0.9530848 | 3_Transplant1 | Podarcis_muralis | Control |
| AD38 | 0.9247925 | 3_Transplant1 | Podarcis_muralis | Control |
| AD39 | 1.0000000 | 3_Transplant1 | Podarcis_muralis | Control |
| AD40 | 0.9991552 | 3_Transplant1 | Podarcis_muralis | Control |
| AD41 | 0.9550964 | 3_Transplant1 | Podarcis_muralis | Control |
| AD42 | 1.0000000 | 3_Transplant1 | Podarcis_muralis | Control |
| AD43 | 1.0000000 | 3_Transplant1 | Podarcis_muralis | Control |
| AD44 | 0.9108807 | 3_Transplant1 | Podarcis_muralis | Control |
| AD46 | 0.8708468 | 3_Transplant1 | Podarcis_muralis | Treatment |
| AD47 | 0.8328653 | 3_Transplant1 | Podarcis_muralis | Treatment |
| AD49 | 0.8940294 | 3_Transplant1 | Podarcis_muralis | Treatment |
| AD50 | 0.8757711 | 3_Transplant1 | Podarcis_muralis | Treatment |
| AD51 | 0.8479193 | 3_Transplant1 | Podarcis_muralis | Treatment |
| AD52 | 0.8596870 | 3_Transplant1 | Podarcis_muralis | Treatment |
| AD53 | 0.8528974 | 3_Transplant1 | Podarcis_muralis | Treatment |
| AD54 | 1.0000000 | 3_Transplant1 | Podarcis_liolepis | Hot_control |
| AD55 | 1.0000000 | 4_Transplant2 | Podarcis_muralis | Control |
| AD56 | 0.9327538 | 4_Transplant2 | Podarcis_muralis | Control |
| AD57 | 1.0000000 | 4_Transplant2 | Podarcis_muralis | Control |
| AD58 | 1.0000000 | 4_Transplant2 | Podarcis_muralis | Control |
| AD59 | 1.0000000 | 4_Transplant2 | Podarcis_muralis | Control |
| AD60 | 0.7988844 | 4_Transplant2 | Podarcis_muralis | Control |
| AD61 | 1.0000000 | 4_Transplant2 | Podarcis_muralis | Control |
| AD62 | 1.0000000 | 4_Transplant2 | Podarcis_muralis | Control |
| AD63 | 0.8992048 | 4_Transplant2 | Podarcis_muralis | Control |
| AD65 | 1.0000000 | 4_Transplant2 | Podarcis_muralis | Treatment |
| AD66 | 0.8821027 | 4_Transplant2 | Podarcis_muralis | Treatment |
| AD68 | 0.8820274 | 4_Transplant2 | Podarcis_muralis | Treatment |
| AD70 | 0.9590793 | 4_Transplant2 | Podarcis_muralis | Treatment |
| AD71 | 0.8360540 | 4_Transplant2 | Podarcis_muralis | Treatment |
| AD72 | 0.8835253 | 4_Transplant2 | Podarcis_muralis | Treatment |
| AD73 | 1.0000000 | 4_Transplant2 | Podarcis_liolepis | Hot_control |
| AD74 | 0.9017271 | 5_Post-FMT1 | Podarcis_muralis | Control |
| AD75 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Treatment |
| AD76 | 0.9862581 | 5_Post-FMT1 | Podarcis_muralis | Treatment |
| AD77 | 0.8963376 | 5_Post-FMT1 | Podarcis_muralis | Treatment |
| AD78 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Treatment |
| AD79 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Control |
| AD80 | 0.8808155 | 5_Post-FMT1 | Podarcis_muralis | Treatment |
| AD81 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Control |
| AD82 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Control |
| AD83 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Treatment |
| AD84 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Treatment |
| AD85 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Control |
| AD86 | 0.9889084 | 5_Post-FMT1 | Podarcis_muralis | Control |
| AD87 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Control |
| AD88 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Treatment |
| AD89 | 0.9562489 | 5_Post-FMT1 | Podarcis_muralis | Control |
| AD90 | 1.0000000 | 5_Post-FMT1 | Podarcis_muralis | Control |
| AD93 | 1.0000000 | 5_Post-FMT1 | Podarcis_liolepis | Hot_control |
| AD94 | 0.7947750 | 5_Post-FMT1 | Podarcis_liolepis | Hot_control |
| AD95 | 0.9506977 | 5_Post-FMT1 | Podarcis_liolepis | Hot_control |
| AD96 | 0.9171413 | 5_Post-FMT1 | Podarcis_liolepis | Hot_control |
| AD97 | 0.9257685 | 5_Post-FMT1 | Podarcis_liolepis | Hot_control |
| AD98 | 0.9048349 | 5_Post-FMT1 | Podarcis_liolepis | Hot_control |
| AD99 | 0.9966152 | 5_Post-FMT1 | Podarcis_liolepis | Hot_control |
| AE01 | 1.0000000 | 5_Post-FMT1 | Podarcis_liolepis | Hot_control |
| AE02 | 1.0000000 | 5_Post-FMT1 | Podarcis_liolepis | Hot_control |
| AE04 | 1.0000000 | 6_Post-FMT2 | Podarcis_muralis | Treatment |
| AE05 | 1.0000000 | 6_Post-FMT2 | Podarcis_muralis | Treatment |
| AE06 | 1.0000000 | 6_Post-FMT2 | Podarcis_muralis | Control |
| AE07 | 0.7797794 | 6_Post-FMT2 | Podarcis_muralis | Control |
| AE08 | 1.0000000 | 6_Post-FMT2 | Podarcis_muralis | Control |
| AE09 | 1.0000000 | 6_Post-FMT2 | Podarcis_muralis | Control |
| AE91 | 0.9259146 | 6_Post-FMT2 | Podarcis_muralis | Treatment |
| AE92 | 1.0000000 | 6_Post-FMT2 | Podarcis_muralis | Control |
| AE93 | 1.0000000 | 6_Post-FMT2 | Podarcis_muralis | Control |
| AE94 | 0.9152960 | 6_Post-FMT2 | Podarcis_muralis | Treatment |
| AE95 | 1.0000000 | 6_Post-FMT2 | Podarcis_muralis | Treatment |
| AE96 | 0.7786112 | 6_Post-FMT2 | Podarcis_muralis | Treatment |
| AE97 | 0.9365037 | 6_Post-FMT2 | Podarcis_muralis | Control |
| AE98 | 1.0000000 | 6_Post-FMT2 | Podarcis_muralis | Control |
| AE99 | 0.8513354 | 6_Post-FMT2 | Podarcis_muralis | Treatment |
| AF01 | 0.9046369 | 6_Post-FMT2 | Podarcis_muralis | Treatment |
| AF02 | 1.0000000 | 6_Post-FMT2 | Podarcis_muralis | Control |
| AF03 | 0.9030087 | 6_Post-FMT2 | Podarcis_muralis | Treatment |
| AF04 | 1.0000000 | 6_Post-FMT2 | Podarcis_liolepis | Hot_control |
| AF05 | 1.0000000 | 6_Post-FMT2 | Podarcis_liolepis | Hot_control |
| AF06 | 1.0000000 | 6_Post-FMT2 | Podarcis_liolepis | Hot_control |
| AF07 | 0.9526443 | 6_Post-FMT2 | Podarcis_liolepis | Hot_control |
| AF08 | 1.0000000 | 6_Post-FMT2 | Podarcis_liolepis | Hot_control |
| AF09 | 1.0000000 | 6_Post-FMT2 | Podarcis_liolepis | Hot_control |
| AF10 | 0.9716162 | 6_Post-FMT2 | Podarcis_liolepis | Hot_control |
| AF11 | 0.9654152 | 6_Post-FMT2 | Podarcis_liolepis | Hot_control |
| AF13 | 0.8818742 | 6_Post-FMT2 | Podarcis_liolepis | Hot_control |
| AFU87 | 0.8810513 | 0_Wild | Podarcis_muralis | Control |
| AFU88 | 0.8317800 | 0_Wild | Podarcis_muralis | Treatment |
| AFU91 | 0.8923699 | 0_Wild | Podarcis_muralis | Treatment |
| AFU92 | 0.8094776 | 0_Wild | Podarcis_muralis | Treatment |
| AFU93 | 0.8517368 | 0_Wild | Podarcis_muralis | Control |
| AFU94 | 0.8325385 | 0_Wild | Podarcis_muralis | Treatment |
| AFU95 | 0.8419270 | 0_Wild | Podarcis_muralis | Treatment |
| AFU96 | 0.8326820 | 0_Wild | Podarcis_muralis | Control |
| AFU97 | 0.8107271 | 0_Wild | Podarcis_muralis | Treatment |
| AFU98 | 0.7506522 | 0_Wild | Podarcis_muralis | Control |
| AFU99 | 0.8582371 | 0_Wild | Podarcis_muralis | Treatment |
| AFV01 | 0.9331539 | 0_Wild | Podarcis_muralis | Treatment |
| AFV02 | 0.8316460 | 0_Wild | Podarcis_muralis | Treatment |
| AFV03 | 0.8752591 | 0_Wild | Podarcis_muralis | Control |
| AFV04 | 0.9180285 | 0_Wild | Podarcis_muralis | Control |
| AFV05 | 1.0000000 | 0_Wild | Podarcis_muralis | Control |
| AFV06 | 1.0000000 | 0_Wild | Podarcis_muralis | Control |
| AFV07 | 0.8460805 | 0_Wild | Podarcis_muralis | Control |
| AFV08 | 0.7497043 | 0_Wild | Podarcis_liolepis | Hot_control |
| AFV09 | 0.5412999 | 0_Wild | Podarcis_liolepis | Hot_control |
| AFV10 | 0.8002499 | 0_Wild | Podarcis_liolepis | Hot_control |
| AFV11 | 0.8225298 | 0_Wild | Podarcis_liolepis | Hot_control |
| AFV12 | 0.7925988 | 0_Wild | Podarcis_liolepis | Hot_control |
| AFV14 | 0.8106269 | 0_Wild | Podarcis_liolepis | Hot_control |
| AFV15 | 0.9691106 | 0_Wild | Podarcis_liolepis | Hot_control |
| AFV16 | 0.8218990 | 0_Wild | Podarcis_liolepis | Hot_control |
| AFV17 | 0.8091152 | 0_Wild | Podarcis_liolepis | Hot_control |
4 MAG catalogue
4.1 Genome phylogeny
# Generate the phylum color heatmap
phylum_heatmap <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
arrange(match(genome, genome_tree$tip.label)) %>%
select(genome,phylum) %>%
mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
column_to_rownames(var = "genome")
# Generate basal tree
circular_tree <- force.ultrametric(genome_tree, method="extend") %>% # extend to ultrametric for the sake of visualisation
ggtree(., layout="fan", open.angle=10, size=0.5)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
# Add phylum ring
circular_tree <- gheatmap(circular_tree, phylum_heatmap, offset=0.55, width=0.1, colnames=FALSE) +
scale_fill_manual(values=phylum_colors) +
geom_tiplab2(size=1, hjust=-0.1) +
theme(legend.position = "none", plot.margin = margin(0, 0, 0, 0), panel.margin = margin(0, 0, 0, 0))
# Flush color scale to enable a new color scheme in the next ring
circular_tree <- circular_tree + new_scale_fill()
# Add completeness ring
circular_tree <- circular_tree +
new_scale_fill() +
scale_fill_gradient(low = "#d1f4ba", high = "#f4baba") +
geom_fruit(
data=genome_metadata,
geom=geom_bar,
mapping = aes(x=completeness, y=genome, fill=contamination),
offset = 0.55,
orientation="y",
stat="identity")
# Add genome-size ring
circular_tree <- circular_tree +
new_scale_fill() +
scale_fill_manual(values = "#cccccc") +
geom_fruit(
data=genome_metadata,
geom=geom_bar,
mapping = aes(x=length, y=genome),
offset = 0.05,
orientation="y",
stat="identity")
# Add text
circular_tree <- circular_tree +
annotate('text', x=3.4, y=0, label=' Phylum', family='ArialMT', size=3.5) +
annotate('text', x=4.2, y=0, label=' Genome quality', family='ArialMT', size=3.5) +
annotate('text', x=4.7, y=0, label=' Genome size', family='ArialMT', size=3.5)
#Plot circular tree
circular_tree %>% open_tree(30) %>% rotate_tree(90)4.2 Genome quality
[1] 92.98658
[1] 7.170467
[1] 2.154888
[1] 2.194372
#Generate quality biplot
genome_biplot <- genome_metadata %>%
select(c(genome,domain,phylum,completeness,contamination,length)) %>%
arrange(match(genome, rev(genome_tree$tip.label))) %>% #sort MAGs according to phylogenetic tree
ggplot(aes(x=completeness,y=contamination,size=length,color=phylum)) +
geom_point(alpha=0.7) +
ylim(c(10,0)) +
scale_color_manual(values=phylum_colors) +
labs(y= "Contamination", x = "Completeness") +
theme_classic() +
theme(legend.position = "none")
#Generate contamination boxplot
genome_contamination <- genome_metadata %>%
ggplot(aes(y=contamination)) +
ylim(c(10,0)) +
geom_boxplot(colour = "#999999", fill="#cccccc") +
theme_void() +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
plot.margin = unit(c(0, 0, 0.40, 0),"inches")) #add bottom-margin (top, right, bottom, left)
#Generate completeness boxplot
genome_completeness <- genome_metadata %>%
ggplot(aes(x=completeness)) +
xlim(c(50,100)) +
geom_boxplot(colour = "#999999", fill="#cccccc") +
theme_void() +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
plot.margin = unit(c(0, 0, 0, 0.50),"inches")) #add left-margin (top, right, bottom, left)
#Render composite figure
grid.arrange(grobs = list(genome_completeness,genome_biplot,genome_contamination),
layout_matrix = rbind(c(1,1,1,1,1,1,1,1,1,1,1,4),
c(2,2,2,2,2,2,2,2,2,2,2,3),
c(2,2,2,2,2,2,2,2,2,2,2,3),
c(2,2,2,2,2,2,2,2,2,2,2,3),
c(2,2,2,2,2,2,2,2,2,2,2,3),
c(2,2,2,2,2,2,2,2,2,2,2,3),
c(2,2,2,2,2,2,2,2,2,2,2,3),
c(2,2,2,2,2,2,2,2,2,2,2,3),
c(2,2,2,2,2,2,2,2,2,2,2,3),
c(2,2,2,2,2,2,2,2,2,2,2,3),
c(2,2,2,2,2,2,2,2,2,2,2,3),
c(2,2,2,2,2,2,2,2,2,2,2,3)))4.3 Functional overview
# Aggregate basal GIFT into elements
function_table <- genome_gifts %>%
to.elements(., GIFT_db)
# Generate basal tree
function_tree <- force.ultrametric(genome_tree, method="extend") %>%
ggtree(., size = 0.3) ***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
#Add phylum colors next to the tree tips
function_tree <- gheatmap(function_tree, phylum_heatmap, offset=0, width=0.1, colnames=FALSE) +
scale_fill_manual(values=phylum_colors) +
labs(fill="Phylum")
#Reset fill scale to use a different colour profile in the heatmap
function_tree <- function_tree + new_scale_fill()
#Add functions heatmap
function_tree <- gheatmap(function_tree, function_table, offset=0.5, width=3.5, colnames=FALSE) +
vexpand(.08) +
coord_cartesian(clip = "off") +
scale_fill_gradient(low = "#f4f4f4", high = "steelblue", na.value="white") +
labs(fill="GIFT")
#Reset fill scale to use a different colour profile in the heatmap
function_tree <- function_tree + new_scale_fill()
# Add completeness barplots
function_tree <- function_tree +
geom_fruit(data=genome_metadata,
geom=geom_bar,
grid.params=list(axis="x", text.size=2, nbreak = 1),
axis.params=list(vline=TRUE),
mapping = aes(x=length, y=genome, fill=completeness),
offset = 3.8,
orientation="y",
stat="identity") +
scale_fill_gradient(low = "#cf8888", high = "#a2cc87") +
labs(fill="Genome\ncompleteness")
function_tree4.4 Functional ordination
# Generate the tSNE ordination
tSNE_function <- Rtsne(X=function_table, dims = 2, check_duplicates = FALSE)
# Plot the ordination
function_ordination <- tSNE_function$Y %>%
as.data.frame() %>%
mutate(genome=rownames(function_table)) %>%
inner_join(genome_metadata, by="genome") %>%
rename(tSNE1="V1", tSNE2="V2") %>%
select(genome,phylum,tSNE1,tSNE2, length) %>%
ggplot(aes(x = tSNE1, y = tSNE2, color = phylum, size=length))+
geom_point(shape=16, alpha=0.7) +
scale_color_manual(values=phylum_colors) +
theme_minimal() +
labs(color="Phylum", size="Genome size") +
guides(color = guide_legend(override.aes = list(size = 5))) # enlarge Phylum dots in legend
function_ordination5 Community composition
5.1 Taxonomy overview
5.1.1 Stacked barplot
# Merge data frames based on sample
transplants_metadata<-sample_metadata%>%
mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)
merged_data<-genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
left_join(., transplants_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
filter(count > 0) #filter 0 counts
ggplot(merged_data, aes(x=newID,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
scale_fill_manual(values=phylum_colors) +
facet_nested(. ~ time_point + type , scales="free") + #facet per day and treatment
guides(fill = guide_legend(ncol = 1)) +
labs(fill="Phylum",y = "Relative abundance",x="Sample")+
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size=0))5.1.1.1 Wild samples
merged_data %>%
filter(time_point=="0_Wild") %>%
ggplot(aes(x=newID,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
scale_fill_manual(values=phylum_colors) +
facet_nested(. ~ Population, scales="free") + #facet per day and treatment
guides(fill = guide_legend(ncol = 1)) +
labs(fill="Phylum",y = "Relative abundance",x="Sample")+
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size=0),
strip.text.x = element_text(size = 12)
)5.1.2 Phylum relative abundances
phylum_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
left_join(genome_metadata, by = join_by(genome == genome)) %>%
group_by(sample,phylum) %>%
summarise(relabun=sum(count))
phylum_summary %>%
group_by(phylum) %>%
summarise(total_mean=mean(relabun*100, na.rm=T),
total_sd=sd(relabun*100, na.rm=T)) %>%
mutate(total=str_c(round(total_mean,2),"±",round(total_sd,2))) %>%
arrange(-total_mean) %>%
dplyr::select(phylum,total) %>%
tt()| phylum | total |
|---|---|
| p__Bacteroidota | 38.75±19.94 |
| p__Bacillota_A | 24.77±15.73 |
| p__Bacillota | 11.96±14.79 |
| p__Pseudomonadota | 9.41±15.86 |
| p__Campylobacterota | 5.49±9.42 |
| p__Verrucomicrobiota | 2.78±6.7 |
| p__Desulfobacterota | 2.37±3.68 |
| p__Chlamydiota | 1.1±6.08 |
| p__Fusobacteriota | 1.06±2.86 |
| p__Cyanobacteriota | 0.93±1.66 |
| p__Bacillota_C | 0.48±0.67 |
| p__Spirochaetota | 0.41±1.25 |
| p__Bacillota_B | 0.26±0.5 |
| p__Actinomycetota | 0.13±0.65 |
| p__Elusimicrobiota | 0.11±0.63 |
5.1.2.1 Wild samples
#Cold and wet population
phylum_summary %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
group_by(phylum, Population) %>%
filter(Population=="Cold_wet") %>%
filter(time_point=="0_Wild") %>%
summarise(total_mean=mean(relabun*100, na.rm=T),
total_sd=sd(relabun*100, na.rm=T)) %>%
mutate(total=str_c(round(total_mean,2),"±",round(total_sd,2))) %>%
arrange(-total_mean) %>%
dplyr::select(phylum,total, Population) %>%
tt()| phylum | total | Population |
|---|---|---|
| p__Bacteroidota | 41.93±15.4 | Cold_wet |
| p__Bacillota_A | 33.96±16.7 | Cold_wet |
| p__Bacillota | 5.73±4.73 | Cold_wet |
| p__Campylobacterota | 5.39±4.64 | Cold_wet |
| p__Verrucomicrobiota | 4.36±3.06 | Cold_wet |
| p__Pseudomonadota | 3.75±3.76 | Cold_wet |
| p__Desulfobacterota | 2.31±1.6 | Cold_wet |
| p__Bacillota_C | 0.85±0.82 | Cold_wet |
| p__Cyanobacteriota | 0.64±0.83 | Cold_wet |
| p__Bacillota_B | 0.5±0.77 | Cold_wet |
| p__Fusobacteriota | 0.37±1.28 | Cold_wet |
| p__Spirochaetota | 0.12±0.35 | Cold_wet |
| p__Actinomycetota | 0.09±0.22 | Cold_wet |
| p__Chlamydiota | 0±0 | Cold_wet |
| p__Elusimicrobiota | 0±0 | Cold_wet |
#Hot and dry population
phylum_summary %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
group_by(phylum, Population) %>%
filter(Population=="Hot_dry")%>%
filter(time_point=="0_Wild") %>%
summarise(total_mean=mean(relabun*100, na.rm=T),
total_sd=sd(relabun*100, na.rm=T)) %>%
mutate(total=str_c(round(total_mean,2),"±",round(total_sd,2))) %>%
arrange(-total_mean) %>%
dplyr::select(phylum,total, Population) %>%
tt()| phylum | total | Population |
|---|---|---|
| p__Bacillota_A | 33.9±16.17 | Hot_dry |
| p__Bacteroidota | 27.63±17.64 | Hot_dry |
| p__Campylobacterota | 13.53±20.98 | Hot_dry |
| p__Pseudomonadota | 10.92±9.94 | Hot_dry |
| p__Bacillota | 6.09±7.66 | Hot_dry |
| p__Spirochaetota | 2.66±2.78 | Hot_dry |
| p__Desulfobacterota | 1.69±1.85 | Hot_dry |
| p__Fusobacteriota | 1.34±1.73 | Hot_dry |
| p__Bacillota_B | 0.76±0.73 | Hot_dry |
| p__Bacillota_C | 0.73±0.5 | Hot_dry |
| p__Cyanobacteriota | 0.52±0.65 | Hot_dry |
| p__Chlamydiota | 0.12±0.18 | Hot_dry |
| p__Verrucomicrobiota | 0.06±0.1 | Hot_dry |
| p__Elusimicrobiota | 0.05±0.14 | Hot_dry |
| p__Actinomycetota | 0±0 | Hot_dry |
phylum_arrange <- phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun)) %>%
arrange(-mean) %>%
select(phylum) %>%
pull()
phylum_summary %>%
filter(phylum %in% phylum_arrange) %>%
mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
geom_jitter(alpha=0.5) +
theme_minimal() +
theme(legend.position="none") +
labs(y="Phylum",x="Relative abundance")5.2 Taxonomy boxplot
5.2.1 Family
family_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,family) %>%
summarise(relabun=sum(count))
family_summary %>%
group_by(family) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| family | mean | sd |
|---|---|---|
| f__Bacteroidaceae | 2.260146e-01 | 0.1384706235 |
| f__Lachnospiraceae | 1.410833e-01 | 0.1062953893 |
| f__Tannerellaceae | 1.045659e-01 | 0.0799894745 |
| f__Helicobacteraceae | 5.448546e-02 | 0.0937279764 |
| f__Mycoplasmoidaceae | 3.756572e-02 | 0.0767893776 |
| f__Erysipelotrichaceae | 3.536287e-02 | 0.0452140267 |
| f__UBA3700 | 3.456595e-02 | 0.0565495010 |
| f__Marinifilaceae | 2.794365e-02 | 0.0272474083 |
| f__Rikenellaceae | 2.725202e-02 | 0.0471735191 |
| f__Enterobacteriaceae | 2.687327e-02 | 0.0929254305 |
| f__Coprobacillaceae | 2.627823e-02 | 0.0907456387 |
| f__ | 2.465083e-02 | 0.0781013121 |
| f__Desulfovibrionaceae | 2.373771e-02 | 0.0367718116 |
| f__DTU072 | 2.183007e-02 | 0.0377159876 |
| f__Ruminococcaceae | 1.832093e-02 | 0.0428115194 |
| f__Rhizobiaceae | 1.579679e-02 | 0.0779688169 |
| f__LL51 | 1.556592e-02 | 0.0616955422 |
| f__UBA3830 | 1.512118e-02 | 0.0441855701 |
| f__Akkermansiaceae | 1.224165e-02 | 0.0317653885 |
| f__Chlamydiaceae | 1.096188e-02 | 0.0607502761 |
| f__Fusobacteriaceae | 1.055726e-02 | 0.0286383947 |
| f__CAG-239 | 9.138683e-03 | 0.0152490206 |
| f__Enterococcaceae | 8.437601e-03 | 0.0473906561 |
| f__Gastranaerophilaceae | 7.848357e-03 | 0.0146292952 |
| f__Oscillospiraceae | 6.624721e-03 | 0.0075288565 |
| f__UBA1997 | 6.613196e-03 | 0.0315103296 |
| f__Streptococcaceae | 6.600789e-03 | 0.0348230465 |
| f__UBA1242 | 4.266475e-03 | 0.0147768596 |
| f__Brevinemataceae | 4.098862e-03 | 0.0125062564 |
| f__Acutalibacteraceae | 3.498766e-03 | 0.0111374416 |
| f__RUG11792 | 2.921450e-03 | 0.0255676374 |
| f__Clostridiaceae | 2.855351e-03 | 0.0174153876 |
| f__UBA660 | 2.600647e-03 | 0.0118148140 |
| f__Peptococcaceae | 2.556218e-03 | 0.0049947786 |
| f__Acidaminococcaceae | 1.980431e-03 | 0.0051045211 |
| f__CAG-508 | 1.874529e-03 | 0.0065256902 |
| f__MGBC116941 | 1.783700e-03 | 0.0077801927 |
| f__Moraxellaceae | 1.540093e-03 | 0.0099192011 |
| f__RUG14156 | 1.428152e-03 | 0.0045670616 |
| f__Staphylococcaceae | 1.411833e-03 | 0.0051727635 |
| f__Anaerovoracaceae | 1.410739e-03 | 0.0027876311 |
| f__Elusimicrobiaceae | 1.088209e-03 | 0.0062780187 |
| f__CAG-288 | 9.840222e-04 | 0.0061275213 |
| f__Anaerotignaceae | 9.320656e-04 | 0.0041174457 |
| f__CALVMC01 | 7.793540e-04 | 0.0044385554 |
| f__Eggerthellaceae | 6.643755e-04 | 0.0021620275 |
| f__Massilibacillaceae | 6.322621e-04 | 0.0016561037 |
| f__Mycobacteriaceae | 6.168531e-04 | 0.0061497354 |
| f__UBA1820 | 4.705627e-04 | 0.0013078764 |
| f__CAG-274 | 4.686117e-04 | 0.0022415212 |
| f__Arcobacteraceae | 3.928587e-04 | 0.0050156837 |
| f__Burkholderiaceae_C | 3.835606e-04 | 0.0048969735 |
| f__Muribaculaceae | 3.508548e-04 | 0.0009525792 |
| f__UBA932 | 3.295199e-04 | 0.0011408058 |
| f__Hepatoplasmataceae | 3.099135e-04 | 0.0039567109 |
| f__Rhodobacteraceae | 3.068016e-04 | 0.0039169801 |
| f__Weeksellaceae | 2.873650e-04 | 0.0032049404 |
| f__Eubacteriaceae | 1.707442e-04 | 0.0006844943 |
| f__Sphingobacteriaceae | 1.561202e-04 | 0.0012685229 |
| f__Devosiaceae | 1.544841e-04 | 0.0015368528 |
| f__Pumilibacteraceae | 1.324439e-04 | 0.0007783049 |
| f__WRAU01 | 9.956857e-05 | 0.0012712064 |
| f__Peptostreptococcaceae | 2.371535e-05 | 0.0003027773 |
family_arrange <- family_summary %>%
group_by(family) %>%
summarise(mean=sum(relabun)) %>%
arrange(-mean) %>%
select(family) %>%
pull()
family_summary %>%
left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
filter(family %in% family_arrange[1:20]) %>%
mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
scale_color_manual(values=phylum_colors[-8]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")5.2.2 Genus
genus_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,genus) %>%
summarise(relabun=sum(count)) %>%
filter(genus != "g__")
genus_summary %>%
group_by(genus) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| genus | mean | sd |
|---|---|---|
| g__Bacteroides | 1.374441e-01 | 0.0923562611 |
| g__Parabacteroides | 9.843371e-02 | 0.0803813676 |
| g__Phocaeicola | 7.109518e-02 | 0.0799972428 |
| g__Helicobacter_J | 3.115738e-02 | 0.0603033642 |
| g__Mycoplasmoides | 3.115302e-02 | 0.0765633496 |
| g__Odoribacter | 2.605667e-02 | 0.0268723184 |
| g__Roseburia | 2.425344e-02 | 0.0567044550 |
| g__NHYM01 | 2.332808e-02 | 0.0810376677 |
| g__Alistipes | 2.224698e-02 | 0.0287419149 |
| g__Coprobacillus | 2.070109e-02 | 0.0894282233 |
| g__Agrobacterium | 1.579679e-02 | 0.0779688169 |
| g__Akkermansia | 1.224165e-02 | 0.0317653885 |
| g__Fusobacterium_A | 1.046073e-02 | 0.0286441117 |
| g__Kineothrix | 9.107908e-03 | 0.0416218616 |
| g__Proteus | 8.976683e-03 | 0.0694135711 |
| g__Dielma | 8.687357e-03 | 0.0090713197 |
| g__CAG-95 | 8.238073e-03 | 0.0207930753 |
| g__JAAYNV01 | 7.265789e-03 | 0.0179169564 |
| g__Desulfovibrio | 7.219276e-03 | 0.0214990147 |
| g__UBA866 | 7.016767e-03 | 0.0295145125 |
| g__Enterococcus | 6.966137e-03 | 0.0463712943 |
| g__Lactococcus | 6.600789e-03 | 0.0348230465 |
| g__Ureaplasma | 6.412700e-03 | 0.0139267552 |
| g__Parabacteroides_B | 6.132159e-03 | 0.0101543965 |
| g__Lacrimispora | 6.028411e-03 | 0.0098068179 |
| g__CALXRO01 | 5.977964e-03 | 0.0313982647 |
| g__Citrobacter | 5.896711e-03 | 0.0340533686 |
| g__NSJ-61 | 5.745781e-03 | 0.0202234473 |
| g__Breznakia | 5.530147e-03 | 0.0240721461 |
| g__Clostridium_AQ | 5.522246e-03 | 0.0123487838 |
| g__Bilophila | 5.044501e-03 | 0.0089558435 |
| g__Hungatella_A | 4.964136e-03 | 0.0096921078 |
| g__Escherichia | 4.342538e-03 | 0.0270859242 |
| g__Salmonella | 4.319018e-03 | 0.0148769561 |
| g__UMGS1251 | 4.312965e-03 | 0.0073601071 |
| g__MGBC136627 | 4.305492e-03 | 0.0164533523 |
| g__Hungatella | 4.150386e-03 | 0.0194068227 |
| g__Clostridium_Q | 4.146767e-03 | 0.0052575243 |
| g__Brevinema | 4.098862e-03 | 0.0125062564 |
| g__Thomasclavelia | 4.046233e-03 | 0.0110779301 |
| g__Scatousia | 3.752075e-03 | 0.0104403539 |
| g__Mailhella | 3.745039e-03 | 0.0104110785 |
| g__Copromonas | 3.643508e-03 | 0.0050495456 |
| g__Enterocloster | 3.613702e-03 | 0.0047492729 |
| g__Ventrimonas | 3.566172e-03 | 0.0071931788 |
| g__Fournierella | 3.313097e-03 | 0.0063192740 |
| g__Limenecus | 3.230504e-03 | 0.0066725343 |
| g__Mucinivorans | 3.006847e-03 | 0.0379999623 |
| g__Lawsonia | 2.916613e-03 | 0.0103686789 |
| g__MGBC133411 | 2.902785e-03 | 0.0074333461 |
| g__Caccovivens | 2.887473e-03 | 0.0112659902 |
| g__Sarcina | 2.855351e-03 | 0.0174153876 |
| g__Eisenbergiella | 2.796704e-03 | 0.0069384489 |
| g__Bacteroides_G | 2.781473e-03 | 0.0352463088 |
| g__CAJLXD01 | 2.730769e-03 | 0.0088951735 |
| g__Acetatifactor | 2.654208e-03 | 0.0055286194 |
| g__Blautia | 2.598789e-03 | 0.0062369300 |
| g__Velocimicrobium | 2.235984e-03 | 0.0067748392 |
| g__C-19 | 2.235603e-03 | 0.0048296119 |
| g__CAZU01 | 2.189719e-03 | 0.0066369837 |
| g__Negativibacillus | 2.145239e-03 | 0.0056002700 |
| g__Intestinimonas | 2.003816e-03 | 0.0035552824 |
| g__Rikenella | 1.998193e-03 | 0.0037323264 |
| g__Phascolarctobacterium | 1.980431e-03 | 0.0051045211 |
| g__Butyricimonas | 1.886974e-03 | 0.0042483569 |
| g__RGIG6463 | 1.855727e-03 | 0.0040258495 |
| g__MGBC116941 | 1.783700e-03 | 0.0077801927 |
| g__JALFVM01 | 1.712574e-03 | 0.0038669765 |
| g__Oscillibacter | 1.546231e-03 | 0.0025273862 |
| g__Acinetobacter | 1.540093e-03 | 0.0099192011 |
| g__Pseudoflavonifractor | 1.489286e-03 | 0.0027026675 |
| g__Citrobacter_A | 1.444197e-03 | 0.0061395639 |
| g__Staphylococcus | 1.411833e-03 | 0.0051727635 |
| g__14-2 | 1.228249e-03 | 0.0098038984 |
| g__RGIG4733 | 1.225656e-03 | 0.0038271024 |
| g__Beduini | 1.217163e-03 | 0.0025500013 |
| g__Scatocola | 1.162463e-03 | 0.0045748144 |
| g__Enterococcus_A | 1.123893e-03 | 0.0100947603 |
| g__UBA1436 | 1.088209e-03 | 0.0062780187 |
| g__Faecisoma | 1.057958e-03 | 0.0056043491 |
| g__RGIG9287 | 9.993630e-04 | 0.0094920364 |
| g__CAG-345 | 9.840222e-04 | 0.0061275213 |
| g__Lachnotalea | 9.593025e-04 | 0.0033990676 |
| g__Blautia_A | 9.546974e-04 | 0.0029560949 |
| g__Ruthenibacterium | 8.602962e-04 | 0.0024818365 |
| g__CAG-269 | 8.276280e-04 | 0.0048017072 |
| g__Marseille-P3106 | 8.230475e-04 | 0.0017874580 |
| g__WRHT01 | 6.666234e-04 | 0.0027445999 |
| g__Eggerthella | 6.643755e-04 | 0.0021620275 |
| g__CHH4-2 | 6.371240e-04 | 0.0020328940 |
| g__Corynebacterium | 6.168531e-04 | 0.0061497354 |
| g__Serratia_A | 6.076344e-04 | 0.0077577570 |
| g__Anaerotruncus | 6.058602e-04 | 0.0016447558 |
| g__RUG14156 | 5.735678e-04 | 0.0021869659 |
| g__RGIG1896 | 5.683407e-04 | 0.0051791669 |
| g__IOR16 | 5.574841e-04 | 0.0016418264 |
| g__Faecimonas | 5.146607e-04 | 0.0054508437 |
| g__CAG-56 | 5.096368e-04 | 0.0016613952 |
| g__MGBC140009 | 4.851579e-04 | 0.0024491799 |
| g__CALURL01 | 4.805911e-04 | 0.0017020401 |
| g__Merdimorpha | 4.705627e-04 | 0.0013078764 |
| g__RGIG8482 | 4.560993e-04 | 0.0030287706 |
| g__Enterobacter | 4.223379e-04 | 0.0042068345 |
| g__Klebsiella | 4.203682e-04 | 0.0049802041 |
| g__Caccenecus | 4.086273e-04 | 0.0018112589 |
| g__Aliarcobacter | 3.928587e-04 | 0.0050156837 |
| g__Scatenecus | 3.851876e-04 | 0.0018282510 |
| g__Alcaligenes | 3.835606e-04 | 0.0048969735 |
| g__Plesiomonas | 3.766988e-04 | 0.0027593254 |
| g__JAHHSE01 | 3.529590e-04 | 0.0014998851 |
| g__HGM05232 | 3.508548e-04 | 0.0009525792 |
| g__Enterococcus_B | 3.475714e-04 | 0.0022665993 |
| g__Egerieousia | 3.295199e-04 | 0.0011408058 |
| g__Stoquefichus | 3.137462e-04 | 0.0020871798 |
| g__Hepatoplasma | 3.099135e-04 | 0.0039567109 |
| g__Paracoccus | 3.068016e-04 | 0.0039169801 |
| g__Moheibacter | 2.873650e-04 | 0.0032049404 |
| g__Scatomorpha | 2.738230e-04 | 0.0010358302 |
| g__Emergencia | 2.601331e-04 | 0.0013298673 |
| g__UBA7185 | 2.523935e-04 | 0.0014817660 |
| g__Eubacterium | 1.707442e-04 | 0.0006844943 |
| g__Sphingobacterium | 1.561202e-04 | 0.0012685229 |
| g__Devosia | 1.544841e-04 | 0.0015368528 |
| g__Anaerosporobacter | 1.507638e-04 | 0.0012978048 |
| g__Caccomorpha | 1.434035e-04 | 0.0010730603 |
| g__UBA2658 | 1.355578e-04 | 0.0007332702 |
| g__Protoclostridium | 1.324439e-04 | 0.0007783049 |
| g__Angelakisella | 1.315171e-04 | 0.0009387427 |
| g__Cetobacterium_A | 9.652924e-05 | 0.0008876688 |
| g__Rahnella | 6.708891e-05 | 0.0008565338 |
| g__Peptostreptococcus | 2.371535e-05 | 0.0003027773 |
genus_arrange <- genus_summary %>%
group_by(genus) %>%
summarise(mean=sum(relabun)) %>%
filter(genus != "g__")%>%
arrange(-mean) %>%
select(genus) %>%
mutate(genus= sub("^g__", "", genus)) %>%
pull()
genus_summary %>%
left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
mutate(genus= sub("^g__", "", genus)) %>%
filter(genus %in% genus_arrange[1:20]) %>%
mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")6 Diversity analysis
6.1 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
functional <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, dist = dist) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(functional = 1) %>%
rownames_to_column(var = "sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample)) %>%
full_join(functional, by = join_by(sample == sample))6.1.1 Wild samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.58) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.2 Acclimation samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.58) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.3 Antibiotics samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="2_Antibiotics") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.58) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.4 Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="3_Transplant1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.7) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.5 Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="4_Transplant2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.7) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.6 Post-Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.7) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.7 Post-Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.7) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.2 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, tree = genome_tree)
beta_q1f <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, dist = dist)6.3 Permanovas
6.3.1 1. Are the wild populations similar?
6.3.1.1 Wild: P.muralis vs P.liolepis
wild <- meta %>%
filter(time_point == "0_Wild")
# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))
wild_nmds <- sample_metadata %>%
filter(time_point == "0_Wild")6.3.1.3 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000012 0.000012 0.0012 999 0.975
Residuals 25 0.257281 0.010291
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.979
Podarcis_muralis 0.97302
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.542719 | 0.2095041 | 6.625717 | 0.001 |
| Residual | 25 | 5.820951 | 0.7904959 | NA | NA |
| Total | 26 | 7.363669 | 1.0000000 | NA | NA |
6.3.1.4 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000048 0.0000476 0.0044 999 0.948
Residuals 25 0.270114 0.0108046
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.946
Podarcis_muralis 0.94763
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.918266 | 0.2608511 | 8.822682 | 0.001 |
| Residual | 25 | 5.435610 | 0.7391489 | NA | NA |
| Total | 26 | 7.353876 | 1.0000000 | NA | NA |
6.3.1.5 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03585 0.035847 2.4912 999 0.135
Residuals 25 0.35973 0.014389
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.133
Podarcis_muralis 0.12705
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.3218613 | 0.2162815 | 6.899207 | 0.001 |
| Residual | 25 | 1.1662981 | 0.7837185 | NA | NA |
| Total | 26 | 1.4881594 | 1.0000000 | NA | NA |
6.3.1.6 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.018367 0.018367 1.5597 999 0.203
Residuals 25 0.294402 0.011776
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.209
Podarcis_muralis 0.22328
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0858578 | 0.172879 | 5.225323 | 0.053 |
| Residual | 25 | 0.4107775 | 0.827121 | NA | NA |
| Total | 26 | 0.4966352 | 1.000000 | NA | NA |
beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))6.3.2 2. Effect of acclimation
accli <- meta %>%
filter(time_point == "1_Acclimation")
# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
accli.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(accli))]
identical(sort(colnames(accli.counts)), sort(as.character(rownames(accli))))
accli_nmds <- sample_metadata %>%
filter(time_point == "1_Acclimation")6.3.2.2 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.11796 0.117959 12.963 999 0.003 **
Residuals 25 0.22748 0.009099
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_wet Hot_dry
Cold_wet 0.002
Hot_dry 0.0013711
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| Population | 1 | 1.639807 | 0.179834 | 5.481634 | 0.001 |
| Residual | 25 | 7.478640 | 0.820166 | NA | NA |
| Total | 26 | 9.118447 | 1.000000 | NA | NA |
6.3.2.3 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.07844 0.078443 5.2384 999 0.032 *
Residuals 25 0.37437 0.014975
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_wet Hot_dry
Cold_wet 0.029
Hot_dry 0.030815
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| Population | 1 | 1.947003 | 0.2306127 | 7.493387 | 0.001 |
| Residual | 25 | 6.495736 | 0.7693873 | NA | NA |
| Total | 26 | 8.442739 | 1.0000000 | NA | NA |
6.3.2.4 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.06739 0.067395 2.9532 999 0.104
Residuals 25 0.57052 0.022821
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_wet Hot_dry
Cold_wet 0.095
Hot_dry 0.098068
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| Population | 1 | 0.2441653 | 0.1224638 | 3.488854 | 0.032 |
| Residual | 25 | 1.7496100 | 0.8775362 | NA | NA |
| Total | 26 | 1.9937754 | 1.0000000 | NA | NA |
6.3.2.5 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.02496 0.024955 0.6729 999 0.454
Residuals 25 0.92714 0.037085
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_wet Hot_dry
Cold_wet 0.448
Hot_dry 0.41979
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| Population | 1 | 0.0279454 | 0.0248037 | 0.6358634 | 0.466 |
| Residual | 25 | 1.0987171 | 0.9751963 | NA | NA |
| Total | 26 | 1.1266624 | 1.0000000 | NA | NA |
beta_q0n_nmds_accli <- beta_div_richness_accli$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli_nmds, by = join_by(sample == Tube_code))
beta_q1n_nmds_accli <- beta_div_neutral_accli$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli_nmds, by = join_by(sample == Tube_code))
beta_q1p_nmds_accli <- beta_div_phylo_accli$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli_nmds, by = join_by(sample == Tube_code))
beta_q1f_nmds_accli <- beta_div_func_accli$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli_nmds, by = join_by(sample == Tube_code))6.3.3 3. Comparison between Wild and Acclimation
accli1 <- meta %>%
filter(time_point == "0_Wild" | time_point == "1_Acclimation")
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
accli1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(accli1))]
identical(sort(colnames(accli1.counts)),sort(as.character(rownames(accli1))))
accli1_nmds <- sample_metadata %>%
filter(time_point == "0_Wild" | time_point == "1_Acclimation")6.3.3.1 Number of samples used
[1] 54
beta_div_richness_accli1<-hillpair(data=accli1.counts, q=0)
beta_div_neutral_accli1<-hillpair(data=accli1.counts, q=1)
beta_div_phylo_accli1<-hillpair(data=accli1.counts, q=1, tree=genome_tree)
beta_div_func_accli1<-hillpair(data=accli1.counts, q=1, dist=dist)6.3.3.1.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.05014 0.050145 6.2252 999 0.011 *
Residuals 52 0.41886 0.008055
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
0_Wild 1_Acclimation
0_Wild 0.014
1_Acclimation 0.015808
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 0.6172653 | 0.0360987 | 3.933397 | 0.001 |
| species | 1 | 2.8279677 | 0.1653842 | 18.020647 | 0.001 |
| individual | 25 | 9.5739861 | 0.5599025 | 2.440331 | 0.001 |
| Residual | 26 | 4.0801621 | 0.2386146 | NA | NA |
| Total | 53 | 17.0993812 | 1.0000000 | NA | NA |
6.3.3.1.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.0199 0.0199035 2.1213 999 0.139
Residuals 52 0.4879 0.0093827
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
0_Wild 1_Acclimation
0_Wild 0.142
1_Acclimation 0.15128
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 0.9050519 | 0.0541893 | 6.651487 | 0.001 |
| species | 1 | 3.3236300 | 0.1989999 | 24.426315 | 0.001 |
| individual | 25 | 8.9352276 | 0.5349902 | 2.626702 | 0.001 |
| Residual | 26 | 3.5377576 | 0.2118206 | NA | NA |
| Total | 53 | 16.7016671 | 1.0000000 | NA | NA |
6.3.3.1.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.01334 0.013340 0.6524 999 0.435
Residuals 52 1.06332 0.020449
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
0_Wild 1_Acclimation
0_Wild 0.437
1_Acclimation 0.42294
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 0.2890434 | 0.0766494 | 7.532050 | 0.001 |
| species | 1 | 0.3508889 | 0.0930498 | 9.143655 | 0.001 |
| individual | 25 | 2.1332925 | 0.5657133 | 2.223620 | 0.001 |
| Residual | 26 | 0.9977533 | 0.2645874 | NA | NA |
| Total | 53 | 3.7709782 | 1.0000000 | NA | NA |
6.3.3.1.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.0123 0.012300 0.4817 999 0.482
Residuals 52 1.3277 0.025533
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
0_Wild 1_Acclimation
0_Wild 0.487
1_Acclimation 0.49073
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 0.0448774 | 0.0269021 | 2.355512 | 0.163 |
| species | 1 | 0.0973005 | 0.0583275 | 5.107077 | 0.034 |
| individual | 25 | 1.0306426 | 0.6178264 | 2.163841 | 0.055 |
| Residual | 26 | 0.4953546 | 0.2969440 | NA | NA |
| Total | 53 | 1.6681751 | 1.0000000 | NA | NA |
beta_richness_nmds_accli1 <- beta_div_richness_accli1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli1_nmds, by = c("sample" = "Tube_code"))
beta_neutral_nmds_accli1 <- beta_div_neutral_accli1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli1_nmds, by = c("sample" = "Tube_code"))
beta_phylo_nmds_accli1 <- beta_div_phylo_accli1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli1_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_accli1 <- beta_div_func_accli1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli1_nmds, by = join_by(sample == Tube_code))6.3.4 4. Do the antibiotics work?
6.3.4.1 Acclimation vs antibiotics
treat <- meta %>%
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))
treat_nmds <- sample_metadata %>%
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")6.3.4.2 Number of samples used
[1] 50
beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)6.3.4.2.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.025318 0.0253178 6.021 999 0.023 *
Residuals 48 0.201837 0.0042049
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.027
2_Antibiotics 0.017817
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 1.888584 | 0.0949462 | 6.361098 | 0.001 |
| species | 1 | 2.117109 | 0.1064350 | 7.130814 | 0.001 |
| individual | 25 | 9.353701 | 0.4702455 | 1.260199 | 0.005 |
| Residual | 22 | 6.531709 | 0.3283734 | NA | NA |
| Total | 49 | 19.891103 | 1.0000000 | NA | NA |
6.3.4.2.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.039587 0.039587 6.8387 999 0.013 *
Residuals 48 0.277854 0.005789
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.015
2_Antibiotics 0.011886
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.024181 | 0.1063620 | 9.051981 | 0.001 |
| species | 1 | 2.853103 | 0.1499183 | 12.758858 | 0.001 |
| individual | 25 | 9.234189 | 0.4852168 | 1.651783 | 0.001 |
| Residual | 22 | 4.919584 | 0.2585029 | NA | NA |
| Total | 49 | 19.031057 | 1.0000000 | NA | NA |
6.3.4.2.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.58372 0.58372 35.413 999 0.001 ***
Residuals 48 0.79119 0.01648
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.001
2_Antibiotics 2.9795e-07
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 1.8065206 | 0.2113909 | 18.636551 | 0.001 |
| species | 1 | 0.7903334 | 0.0924813 | 8.153292 | 0.001 |
| individual | 25 | 3.8164689 | 0.4465860 | 1.574869 | 0.005 |
| Residual | 22 | 2.1325541 | 0.2495419 | NA | NA |
| Total | 49 | 8.5458771 | 1.0000000 | NA | NA |
6.3.4.2.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.18591 0.185914 5.0679 999 0.023 *
Residuals 48 1.76088 0.036685
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.032
2_Antibiotics 0.028989
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 1.8020952 | 0.3750193 | 33.6195614 | 0.001 |
| species | 1 | 0.0031247 | 0.0006503 | 0.0582938 | 0.848 |
| individual | 25 | 1.8208629 | 0.3789249 | 1.3587875 | 0.222 |
| Residual | 22 | 1.1792568 | 0.2454055 | NA | NA |
| Total | 49 | 4.8053396 | 1.0000000 | NA | NA |
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_treat <- beta_div_func_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))6.3.5 5. Does the FMT work?
6.3.5.1 Comparison between FMT2 vs Post-FMT2
#Create newID to identify duplicated samples
transplants_metadata<-sample_metadata%>%
mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)
transplant3<-transplants_metadata%>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")%>%
column_to_rownames("newID")
transplant3_nmds <- transplants_metadata %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")
full_counts<-temp_genome_counts %>%
t()%>%
as.data.frame()%>%
rownames_to_column("Tube_code")%>%
full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))
transplant3_counts<-full_counts %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2") %>%
subset(select=-c(315:324)) %>%
column_to_rownames("newID")%>%
subset(select=-c(1))%>%
t() %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric)
identical(sort(colnames(transplant3_counts)),sort(as.character(rownames(transplant3))))6.3.5.2 Number of samples used
[1] 49
beta_div_richness_transplant3<-hillpair(data=transplant3_counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3_counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3_counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3_counts, q=1, dist=dist)
#Arrange of metadata dataframe
transplant3_arrange<-transplant3[labels(beta_div_neutral_transplant3$S),]6.3.5.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.180473 | 0.0855095 | 6.984555 | 0.001 |
| time_point | 1 | 0.860906 | 0.0623612 | 5.093759 | 0.001 |
| type | 1 | 1.459433 | 0.1057165 | 8.635089 | 0.001 |
| individual | 24 | 6.755100 | 0.4893170 | 1.665341 | 0.001 |
| Residual | 21 | 3.549250 | 0.2570959 | NA | NA |
| Total | 48 | 13.805162 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 1.4169018 5.739828 0.15622903 0.001 0.003 *
2 Control vs Hot_control 1 2.0940966 8.509112 0.21005427 0.001 0.003 *
3 Treatment vs Hot_control 1 0.3004618 1.265034 0.04179854 0.159 0.477
6.3.5.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.2800927 | 0.0939787 | 8.796453 | 0.001 |
| time_point | 1 | 0.9350566 | 0.0686477 | 6.425458 | 0.001 |
| type | 1 | 1.9135997 | 0.1404879 | 13.149743 | 0.001 |
| individual | 24 | 6.4363516 | 0.4725281 | 1.842870 | 0.001 |
| Residual | 21 | 3.0559984 | 0.2243577 | NA | NA |
| Total | 48 | 13.6210990 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 1.8758788 8.282671 0.21084796 0.001 0.003 *
2 Control vs Hot_control 1 2.4396317 10.635546 0.24945256 0.001 0.003 *
3 Treatment vs Hot_control 1 0.3158428 1.394345 0.04587515 0.126 0.378
6.3.5.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.1400466 | 0.0952654 | 6.956436 | 0.002 |
| time_point | 1 | 0.1138047 | 0.0774145 | 5.652939 | 0.001 |
| type | 1 | 0.1432667 | 0.0974558 | 7.116383 | 0.001 |
| individual | 24 | 0.6501795 | 0.4422784 | 1.345663 | 0.041 |
| Residual | 21 | 0.4227709 | 0.2875859 | NA | NA |
| Total | 48 | 1.4700683 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.14387705 5.735321 0.15612552 0.001 0.003 *
2 Control vs Hot_control 1 0.22715701 9.044894 0.22036587 0.001 0.003 *
3 Treatment vs Hot_control 1 0.04648319 1.704277 0.05550617 0.129 0.387
6.3.5.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0092808 | 0.0077189 | 0.4182529 | 0.493 |
| time_point | 1 | -0.0061674 | -0.0051295 | -0.2779456 | 0.895 |
| type | 1 | 0.0831052 | 0.0691191 | 3.7452726 | 0.093 |
| individual | 24 | 0.6501528 | 0.5407359 | 1.2208414 | 0.345 |
| Residual | 21 | 0.4659767 | 0.3875556 | NA | NA |
| Total | 48 | 1.2023481 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.078539743 4.59293783 0.129040706 0.069 0.207
2 Control vs Hot_control 1 0.052468954 2.13675422 0.062593948 0.165 0.495
3 Treatment vs Hot_control 1 -0.002340352 -0.07432315 -0.002569452 0.887 1.000
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == newID))
beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == newID))
beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == newID))
beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == newID))p0<-beta_richness_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
p1<-beta_neutral_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylo_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_func_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")6.3.5.3 Comparison between the different experimental time points
6.3.5.3.1 Comparison against Wild samples
The estimated time for calculating the 2850 pairwise combinations is 14 seconds.
6.3.6 6. Are there differences between the control and the treatment group?
6.3.6.2 Number of samples used
[1] 26
beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post1_arrange<-post1[labels(beta_div_neutral_post1$S),]6.3.6.2.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.017675 0.0088373 2.3825 999 0.099 .
Residuals 23 0.085312 0.0037092
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.0060000 0.661
Hot_control 0.0068795 0.212
Treatment 0.6248469 0.2084296
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.6340254 | 0.0768024 | 2.065607 | 0.004 |
| type | 1 | 0.5615418 | 0.0680222 | 1.829461 | 0.010 |
| Residual | 23 | 7.0597099 | 0.8551754 | NA | NA |
| Total | 25 | 8.2552771 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.5615418 1.729004 0.1033537 0.016 0.048 .
2 Control vs Hot_control 1 0.8438429 2.793772 0.1486541 0.002 0.006 *
3 Treatment vs Hot_control 1 0.3734921 1.268929 0.0779971 0.098 0.294
6.3.6.2.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.011001 0.0055005 0.6303 999 0.564
Residuals 23 0.200714 0.0087267
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.21500 0.957
Hot_control 0.21166 0.467
Treatment 0.95468 0.43604
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.7907904 | 0.1076445 | 3.056657 | 0.001 |
| type | 1 | 0.6051778 | 0.0823784 | 2.339205 | 0.009 |
| Residual | 23 | 5.9503501 | 0.8099772 | NA | NA |
| Total | 25 | 7.3463184 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.6051778 2.250849 0.13047758 0.009 0.027 .
2 Control vs Hot_control 1 1.0528902 4.143637 0.20570451 0.001 0.003 *
3 Treatment vs Hot_control 1 0.4150076 1.637268 0.09840968 0.046 0.138
6.3.6.2.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.00440 0.0021994 0.1369 999 0.91
Residuals 23 0.36941 0.0160614
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.92400 0.694
Hot_control 0.91505 0.803
Treatment 0.63312 0.73046
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0560850 | 0.0531376 | 1.3149967 | 0.280 |
| type | 1 | 0.0184254 | 0.0174571 | 0.4320099 | 0.791 |
| Residual | 23 | 0.9809570 | 0.9294053 | NA | NA |
| Total | 25 | 1.0554673 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.01842535 0.4144162 0.02688498 0.775 1.000
2 Control vs Hot_control 1 0.05987967 1.7387847 0.09802164 0.117 0.351
3 Treatment vs Hot_control 1 0.03212966 0.6477782 0.04139746 0.692 1.000
6.3.6.2.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.00400 0.0020014 0.145 999 0.864
Residuals 23 0.31753 0.0138057
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.59900 0.750
Hot_control 0.59817 0.849
Treatment 0.75141 0.83718
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0024979 | 0.0033024 | 0.0900845 | 0.638 |
| type | 1 | 0.1161466 | 0.1535542 | 4.1887855 | 0.063 |
| Residual | 23 | 0.6377435 | 0.8431434 | NA | NA |
| Total | 25 | 0.7563879 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.11614656 4.724791 0.23953568 0.069 0.207
2 Control vs Hot_control 1 0.05000930 1.704826 0.09629160 0.224 0.672
3 Treatment vs Hot_control 1 0.01235859 0.423812 0.02747777 0.503 1.000
beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
p1<-beta_neutral_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")6.3.6.4 Number of samples used
[1] 27
beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post2_arrange<-post2[labels(beta_div_neutral_post2$S),]6.3.6.4.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.002011 0.0010056 0.1982 999 0.83
Residuals 24 0.121775 0.0050740
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.70300 0.786
Hot_control 0.67789 0.624
Treatment 0.79246 0.59820
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | 1.504341 | 0.1967776 | 2.939822 | 0.001 |
| Residual | 24 | 6.140538 | 0.8032224 | NA | NA |
| Total | 26 | 7.644879 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Treatment vs Control 1 0.6463814 2.560441 0.1379515 0.001 0.003 *
2 Treatment vs Hot_control 1 0.4796256 1.916520 0.1069694 0.001 0.003 *
3 Control vs Hot_control 1 1.1305044 4.268317 0.2105906 0.001 0.003 *
6.3.6.4.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.008262 0.0041311 0.8024 999 0.468
Residuals 24 0.123559 0.0051483
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.46700 0.694
Hot_control 0.44675 0.241
Treatment 0.65989 0.25095
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | 1.923807 | 0.2603795 | 4.224537 | 0.001 |
| Residual | 24 | 5.464666 | 0.7396205 | NA | NA |
| Total | 26 | 7.388473 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Treatment vs Control 1 1.0227481 4.648335 0.2251191 0.001 0.003 *
2 Treatment vs Hot_control 1 0.5010202 2.206532 0.1211945 0.001 0.003 *
3 Control vs Hot_control 1 1.3619424 5.771031 0.2650785 0.001 0.003 *
6.3.6.4.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.000407 0.0002034 0.0487 999 0.959
Residuals 24 0.100305 0.0041794
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.92800 0.874
Hot_control 0.93765 0.769
Treatment 0.83933 0.76015
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | 0.1594363 | 0.2042241 | 3.079623 | 0.001 |
| Residual | 24 | 0.6212564 | 0.7957759 | NA | NA |
| Total | 26 | 0.7806927 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Treatment vs Control 1 0.05927454 2.382025 0.1295845 0.021 0.063
2 Treatment vs Hot_control 1 0.06906280 2.722460 0.1454115 0.005 0.015 .
3 Control vs Hot_control 1 0.11081709 4.043656 0.2017424 0.001 0.003 *
6.3.6.4.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.01126 0.0056302 0.2861 999 0.773
Residuals 24 0.47233 0.0196806
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.54800 0.634
Hot_control 0.48255 0.790
Treatment 0.60116 0.75643
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | -0.0038724 | -0.0056213 | -0.0670788 | 0.936 |
| Residual | 24 | 0.6927468 | 1.0056213 | NA | NA |
| Total | 26 | 0.6888744 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Treatment vs Control 1 -0.008527330 -0.46290555 -0.029793572 0.853 1
2 Treatment vs Hot_control 1 -0.001648721 -0.04717131 -0.002956924 0.910 1
3 Control vs Hot_control 1 0.004367477 0.13147026 0.008149924 0.690 1
beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
p1<-beta_neutral_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")7 Functional differences
7.1 Data preparation
# Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts, GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts$genome, ]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
elements <- GIFTs_elements_filtered %>%
as.data.frame()
# Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered, GIFT_db)
functions <- GIFTs_functions %>%
as.data.frame()
# Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions, GIFT_db)
domains <- GIFTs_domains %>%
as.data.frame()
# Get community-weighed average GIFTs per sample
GIFTs_elements_community <- to.community(GIFTs_elements_filtered, genome_counts_filt %>% column_to_rownames(., "genome") %>% tss(), GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions, genome_counts_filt %>% column_to_rownames(., "genome") %>% tss(), GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains, genome_counts_filt %>% column_to_rownames(., "genome") %>% tss(), GIFT_db)
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)7.2 Genomes GIFT profiles
GIFTs_elements %>%
as_tibble(., rownames = "MAG") %>%
reshape2::melt() %>%
rename(Code_element = variable, GIFT = value) %>%
inner_join(GIFT_db,by="Code_element") %>%
ggplot(., aes(x=Code_element, y=MAG, fill=GIFT, group=Code_function))+
geom_tile()+
scale_y_discrete(guide = guide_axis(check.overlap = TRUE))+
scale_x_discrete(guide = guide_axis(check.overlap = TRUE))+
scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
facet_grid(. ~ Code_function, scales = "free", space = "free")+
theme_grey(base_size=8)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text.x = element_text(angle = 90))7.3 Function level
GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(var="sample") %>%
filter(sample!="AD69") %>%
pivot_longer(!sample,names_to="trait",values_to="gift") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
ggplot(aes(x=trait,y=time_point,fill=gift)) +
geom_tile(colour="white", size=0.2)+
scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
facet_grid(type ~ ., scales="free",space="free")7.4 Element level
GIFTs_elements_community_merged<-GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(var="sample") %>%
filter(sample!="AD69") %>%
pivot_longer(!sample,names_to="trait",values_to="gift") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code))%>%
mutate(functionid = substr(trait, 1, 3)) %>%
mutate(trait = case_when(
trait %in% GIFT_db$Code_element ~ GIFT_db$Element[match(trait, GIFT_db$Code_element)],
TRUE ~ trait
)) %>%
mutate(functionid = case_when(
functionid %in% GIFT_db$Code_function ~ GIFT_db$Function[match(functionid, GIFT_db$Code_function)],
TRUE ~ functionid
)) %>%
mutate(trait=factor(trait,levels=unique(GIFT_db$Element))) %>%
mutate(functionid=factor(functionid,levels=unique(GIFT_db$Function)))
# Create an interaction variable for time_point and sample
GIFTs_elements_community_merged$interaction_var <- interaction(GIFTs_elements_community_merged$sample, GIFTs_elements_community_merged$time_point)
ggplot(GIFTs_elements_community_merged,aes(x=interaction_var,y=trait,fill=gift)) +
geom_tile(colour="white", linewidth=0.2)+
scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
facet_grid(functionid ~ type, scales="free",space="free") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=5),
strip.text.y = element_text(angle = 0)) +
labs(y="Traits",x="Time_point",fill="GIFT")+
scale_x_discrete(labels = function(x) gsub(".*\\.", "", x))7.5 Comparison of samples from the 0 Time_point (0_Wild)
7.5.1 GIFTs Functional community
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
group_by(species) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
species MCI sd
<chr> <dbl> <dbl>
1 Podarcis_liolepis 0.327 0.0244
2 Podarcis_muralis 0.346 0.0194
7.5.1.1 GIFT test visualisation
GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
select(c(1:21, 27)) %>%
pivot_longer(-c(sample,type),names_to = "trait", values_to = "value") %>%
mutate(trait = case_when(
trait %in% GIFT_db$Code_function ~ GIFT_db$Function[match(trait, GIFT_db$Code_function)],
TRUE ~ trait
)) %>%
mutate(trait=factor(trait,levels=unique(GIFT_db$Function))) %>%
ggplot(aes(x=value, y=type, group=type, fill=type, color=type)) +
geom_boxplot() +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_grid(trait ~ ., space="free", scales="free") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
strip.text.y = element_text(angle = 0)) +
labs(y="Traits",x="Metabolic capacity index")7.5.2 GIFTs Domain community
GIFTs_domains_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
group_by(species) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
species MCI sd
<chr> <dbl> <dbl>
1 Podarcis_liolepis 0.375 0.0315
2 Podarcis_muralis 0.390 0.0208
7.5.3 GIFTs Elements community
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
group_by(species) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
species MCI sd
<chr> <dbl> <dbl>
1 Podarcis_liolepis 0.313 0.0329
2 Podarcis_muralis 0.345 0.0233
sample_metadata_wild <- sample_metadata%>%
filter(time_point == "0_Wild")
element_gift_wild <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "Tube_code") %>%
inner_join(., sample_metadata_wild[c(1,3)], by="Tube_code")# Find numeric columns
numeric_cols <- sapply(element_gift_wild, is.numeric)
# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_wild[, numeric_cols])
# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]
# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_wild[, !numeric_cols | colnames(element_gift_wild) %in% nonzero_numeric_cols]significant_elements_wild <- filtered_data %>%
pivot_longer(-c(Tube_code,species), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ species, exact=FALSE)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05) %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift_wild %>%
dplyr::select(-c(species)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_wild$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Tube_code")%>%
left_join(., sample_metadata_wild[c(1,3)], by = join_by(Tube_code == Tube_code))
element_gift_filt %>%
dplyr::select(-Tube_code)%>%
group_by(species) %>%
summarise(across(everything(), mean))%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))
element_gift_names <- element_gift_filt%>%
dplyr::select(-species)%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
dplyr::select(-Elements)%>%
dplyr::select(Function, everything())%>%
t()%>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Tube_code")%>%
left_join(., sample_metadata_wild[c(1,4)], by = join_by(Tube_code == Tube_code))colNames <- names(element_gift_names)[2:30] #always check names(element_gift_names) first to know where your traits finish
for(i in colNames){
plt <- ggplot(element_gift_names, aes(x=Population, y=.data[[i]], color = Population)) +
geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
geom_jitter(width = 0.1, show.legend = TRUE) +
theme_minimal() +
theme(
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank())
print(plt)
}7.6 Comparison of samples from the 1st Time_point (1_Acclimation)
7.6.1 GIFTs Functional community
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
group_by(species) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
species MCI sd
<chr> <dbl> <dbl>
1 Podarcis_liolepis 0.348 0.0158
2 Podarcis_muralis 0.331 0.0321
7.6.1.1 GIFT test visualisation
GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
select(c(1:21, 27)) %>%
pivot_longer(-c(sample,type),names_to = "trait", values_to = "value") %>%
mutate(trait = case_when(
trait %in% GIFT_db$Code_function ~ GIFT_db$Function[match(trait, GIFT_db$Code_function)],
TRUE ~ trait
)) %>%
mutate(trait=factor(trait,levels=unique(GIFT_db$Function))) %>%
ggplot(aes(x=value, y=type, group=type, fill=type, color=type)) +
geom_boxplot() +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_grid(trait ~ ., space="free", scales="free") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
strip.text.y = element_text(angle = 0)) +
labs(y="Traits",x="Metabolic capacity index")7.6.2 GIFTs Domain community
GIFTs_domains_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
group_by(species) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
species MCI sd
<chr> <dbl> <dbl>
1 Podarcis_liolepis 0.395 0.0211
2 Podarcis_muralis 0.370 0.0307
7.6.3 GIFTs Elements community
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
group_by(species) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
species MCI sd
<chr> <dbl> <dbl>
1 Podarcis_liolepis 0.350 0.0225
2 Podarcis_muralis 0.332 0.0316
sample_metadata_accli <- sample_metadata%>%
filter(time_point == "1_Acclimation")
element_gift_accli <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "Tube_code") %>%
inner_join(., sample_metadata_accli[c(1,3)], by="Tube_code")# Find numeric columns
numeric_cols <- sapply(element_gift_accli, is.numeric)
# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_accli[, numeric_cols])
# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]
# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_accli[, !numeric_cols | colnames(element_gift_accli) %in% nonzero_numeric_cols]significant_elements_accli <- filtered_data %>%
pivot_longer(-c(Tube_code,species), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ species, exact=FALSE)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05) %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift_accli %>%
dplyr::select(-c(species)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_accli$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Tube_code")%>%
left_join(., sample_metadata_accli[c(1,3)], by = join_by(Tube_code == Tube_code))
element_gift_filt %>%
dplyr::select(-Tube_code)%>%
group_by(species) %>%
summarise(across(everything(), mean))%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))
element_gift_names <- element_gift_filt%>%
dplyr::select(-species)%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
dplyr::select(-Elements)%>%
dplyr::select(Function, everything())%>%
t()%>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Tube_code")%>%
left_join(., sample_metadata_accli[c(1,4)], by = join_by(Tube_code == Tube_code))colNames <- names(element_gift_names)[2:14] #always check names(element_gift_names) first to know where your traits finish
for(i in colNames){
plt <- ggplot(element_gift_names, aes(x=Population, y=.data[[i]], color = Population)) +
geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
geom_jitter(width = 0.1, show.legend = TRUE) +
theme_minimal() +
theme(
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank())
print(plt)
}7.7 Comparison of samples from the 5th Time_point (5_Post-FMT1)
7.7.1 GIFTs Functional community
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
group_by(type) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 3 × 3
type MCI sd
<chr> <dbl> <dbl>
1 Control 0.373 0.0247
2 Hot_control 0.372 0.0367
3 Treatment 0.353 0.0186
7.7.1.1 GIFT test visualisation
GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
select(c(1:21, 27)) %>%
pivot_longer(-c(sample,type),names_to = "trait", values_to = "value") %>%
mutate(trait = case_when(
trait %in% GIFT_db$Code_function ~ GIFT_db$Function[match(trait, GIFT_db$Code_function)],
TRUE ~ trait
)) %>%
mutate(trait=factor(trait,levels=unique(GIFT_db$Function))) %>%
ggplot(aes(x=value, y=type, group=type, fill=type, color=type)) +
geom_boxplot() +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_grid(trait ~ ., space="free", scales="free") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
strip.text.y = element_text(angle = 0)) +
labs(y="Traits",x="Metabolic capacity index")7.7.2 GIFTs Domain community
GIFTs_domains_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
group_by(type) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 3 × 3
type MCI sd
<chr> <dbl> <dbl>
1 Control 0.412 0.0213
2 Hot_control 0.415 0.0437
3 Treatment 0.391 0.0263
7.7.3 GIFTs Elements community
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
group_by(type) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 3 × 3
type MCI sd
<chr> <dbl> <dbl>
1 Control 0.380 0.0280
2 Hot_control 0.379 0.0372
3 Treatment 0.359 0.0214
sample_metadata_tm5 <- sample_metadata%>%
filter(time_point == "5_Post-FMT1")%>%
filter(type != "Hot_control")
element_gift_tm5 <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "Tube_code") %>%
inner_join(sample_metadata_tm5 %>% select(1, 7), by = "Tube_code")# Find numeric columns
numeric_cols <- sapply(element_gift_tm5, is.numeric)
# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_tm5[, numeric_cols])
# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]
# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_tm5[, !numeric_cols | colnames(element_gift_tm5) %in% nonzero_numeric_cols]significant_elements_tm5 <- filtered_data %>%
pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ type, exact=FALSE)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_value < 0.05) %>% #take into account that p_value is used and not p_adjust
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift_tm5 %>%
dplyr::select(-c(type)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_tm5$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Tube_code")%>%
left_join(., sample_metadata_tm5[c(1,7)], by = join_by(Tube_code == Tube_code))
element_gift_filt %>%
dplyr::select(-Tube_code)%>%
group_by(type) %>%
summarise(across(everything(), mean))%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))
element_gift_names <- element_gift_filt%>%
dplyr::select(-type)%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
dplyr::select(-Elements)%>%
dplyr::select(Function, everything())%>%
t()%>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Tube_code")%>%
left_join(., sample_metadata_tm5[c(1,7)], by = join_by(Tube_code == Tube_code))colNames <- names(element_gift_names)[2:14] #always check names(element_gift_names) first to now where your traits finish
for(i in colNames){
plt <- ggplot(element_gift_names, aes(x=type, y=.data[[i]], color = type)) +
geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
geom_jitter(width = 0.1, show.legend = TRUE) +
theme_minimal() +
theme(
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank())
print(plt)
}7.8 Comparison of samples from the 6th Time_point (6_Post-FMT2)
7.8.1 GIFTs Functional community
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
group_by(type) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 3 × 3
type MCI sd
<chr> <dbl> <dbl>
1 Control 0.352 0.0223
2 Hot_control 0.350 0.0293
3 Treatment 0.346 0.0255
7.8.1.1 GIFT test visualisation
GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
select(c(1:21, 27)) %>%
pivot_longer(-c(sample,type),names_to = "trait", values_to = "value") %>%
mutate(trait = case_when(
trait %in% GIFT_db$Code_function ~ GIFT_db$Function[match(trait, GIFT_db$Code_function)],
TRUE ~ trait
)) %>%
mutate(trait=factor(trait,levels=unique(GIFT_db$Function))) %>%
ggplot(aes(x=value, y=type, group=type, fill=type, color=type)) +
geom_boxplot() +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_grid(trait ~ ., space="free", scales="free") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
strip.text.y = element_text(angle = 0)) +
labs(y="Traits",x="Metabolic capacity index")7.8.2 GIFTs Domain community
GIFTs_domains_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
group_by(type) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 3 × 3
type MCI sd
<chr> <dbl> <dbl>
1 Control 0.399 0.0171
2 Hot_control 0.388 0.0321
3 Treatment 0.392 0.0240
7.8.3 GIFTs Elements community
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
group_by(type) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 3 × 3
type MCI sd
<chr> <dbl> <dbl>
1 Control 0.357 0.0215
2 Hot_control 0.347 0.0302
3 Treatment 0.350 0.0293
sample_metadata_TM6 <- sample_metadata%>%
filter(time_point == "6_Post-FMT2")%>%
filter(type != "Hot_control")
element_gift_TM6 <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "Tube_code") %>%
inner_join(sample_metadata_TM6 %>% select(1, 7), by = "Tube_code")# Find numeric columns
numeric_cols <- sapply(element_gift_TM6, is.numeric)
# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_TM6[, numeric_cols])
# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]
# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_TM6[, !numeric_cols | colnames(element_gift_TM6) %in% nonzero_numeric_cols]significant_elements_TM6 <- filtered_data %>%
pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ type, exact=FALSE)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_value < 0.05) %>% #take into account that p_value is used and not p_adjust
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift_TM6 %>%
dplyr::select(-c(type)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_TM6$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Tube_code")%>%
left_join(., sample_metadata_TM6[c(1,7)], by = join_by(Tube_code == Tube_code))
element_gift_filt %>%
dplyr::select(-Tube_code)%>%
group_by(type) %>%
summarise(across(everything(), mean))%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))
element_gift_names <- element_gift_filt%>%
dplyr::select(-type)%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
dplyr::select(-Elements)%>%
dplyr::select(Function, everything())%>%
t()%>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Tube_code")%>%
left_join(., sample_metadata_TM6[c(1,7)], by = join_by(Tube_code == Tube_code))colNames <- names(element_gift_names)[2:20] #always check names(element_gift_names) first to now where your traits finish
for(i in colNames){
plt <- ggplot(element_gift_names, aes(x=type, y=.data[[i]], color = type)) +
geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
geom_jitter(width = 0.1, show.legend = TRUE) +
theme_minimal() +
theme(
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank())
print(plt)
}7.9 Domain level
7.9.1 Comparison of samples from the 0 Time_point (0_Wild)
#Merge the functional domains with the metadata
merge_gift_wild<- GIFTs_domains_community %>%
as.data.frame() %>%
rownames_to_column(., "Tube_code") %>%
inner_join(., sample_metadata_wild, by="Tube_code")#Biosynthesis
p1 <-merge_gift_wild %>%
ggplot(aes(x=species,y=Biosynthesis,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")
#Degradation
p2 <-merge_gift_wild %>%
ggplot(aes(x=species,y=Degradation,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.45, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")
#Structure
p3 <-merge_gift_wild %>%
ggplot(aes(x=species,y=Structure,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 3, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")
#Overall
p4 <-merge_gift_wild %>%
ggplot(aes(x=species,y=Overall,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")7.9.2 Comparison of samples from the 1st Time_point (1_Acclimation)
#Merge the functional domains with the metadata
merge_gift_accli<- GIFTs_domains_community %>%
as.data.frame() %>%
rownames_to_column(., "Tube_code") %>%
inner_join(., sample_metadata_accli, by="Tube_code")#Biosynthesis
p1 <-merge_gift_accli %>%
ggplot(aes(x=species,y=Biosynthesis,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")
#Degradation
p2 <-merge_gift_accli %>%
ggplot(aes(x=species,y=Degradation,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.45, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")
#Structure
p3 <-merge_gift_accli %>%
ggplot(aes(x=species,y=Structure,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 3, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")
#Overall
p4 <-merge_gift_accli %>%
ggplot(aes(x=species,y=Overall,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")7.9.3 Comparison of samples from the 5th Time_point (5_Post-FMT1)
#Merge the functional domains with the metadata
merge_gift_tm5<- GIFTs_domains_community %>%
as.data.frame() %>%
rownames_to_column(., "Tube_code") %>%
inner_join(., sample_metadata_tm5, by="Tube_code")#Biosynthesis
p1 <-merge_gift_tm5 %>%
ggplot(aes(x=species,y=Biosynthesis,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")
#Degradation
p2 <-merge_gift_tm5 %>%
ggplot(aes(x=species,y=Degradation,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.45, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")
#Structure
p3 <-merge_gift_tm5 %>%
ggplot(aes(x=species,y=Structure,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 3, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")
#Overall
p4 <-merge_gift_tm5 %>%
ggplot(aes(x=species,y=Overall,color=species,fill=species))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Species")Warning: Computation failed in `stat_compare_means()`.
Computation failed in `stat_compare_means()`.
Computation failed in `stat_compare_means()`.
Computation failed in `stat_compare_means()`.
Caused by error:
! argument "x" is missing, with no default
7.9.4 Comparison of samples from the 6th Time_point (6_Post-FMT2)
#Merge the functional domains with the metadata
merge_gift_TM6 <- GIFTs_domains_community %>%
as.data.frame() %>%
rownames_to_column(., "Tube_code") %>%
inner_join(., sample_metadata_TM6, by="Tube_code")#Biosynthesis
p1 <-merge_gift_TM6 %>%
ggplot(aes(x=type,y=Biosynthesis,color=type,fill=type))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "Type")
#Degradation
p2 <-merge_gift_TM6 %>%
ggplot(aes(x=type,y=Degradation,color=type,fill=type))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "type")
#Structure
p3 <-merge_gift_TM6 %>%
ggplot(aes(x=type,y=Structure,color=type,fill=type))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 3, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "type")
#Overall
p4 <-merge_gift_TM6 %>%
ggplot(aes(x=type,y=Overall,color=type,fill=type))+
geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+
geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
stat_compare_means() +
theme(axis.text.x = element_text(vjust = 0.5, size=10),
axis.text.y = element_text(size=10),
axis.title=element_text(size=12,face="bold"),
axis.text = element_text(face="bold", size=18),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size=10),
legend.title = element_text(size=12),
legend.position="none",
legend.key.size = unit(1, 'cm'),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
labs( x = "type")University of Copenhagen, garazi.bideguren@sund.ku.dk↩︎
University of Copenhagen, ostaizka.aizpurua@sund.ku.dk↩︎
Sociedad de Ciencias Aranzadi-Departamento de Herpetología, ccabido@aranzadi.eus↩︎
University of Copenhagen, antton.alberdi@sund.ku.dk↩︎